# clear environment
rm(list=ls())

# load packages
library(openxlsx)
library(readstata13)
library(marginaleffects)
library(ggplot2)

# load data set with chapter lengths and crime numbers
dat = openxlsx::read.xlsx("vsberichte_metadata.xlsx", sheet = 1)

# binary varibable indicating whether interior minister is from a center-right party
dat$cr_intmin <- ifelse(dat$intmin_party == "cdu" | dat$intmin_party== "csu" | dat$intmin_party=="fdp", 1, ifelse(dat$intmin_party=="pro", NA, 0))

# remove those without center-right or center-left interior minister
dat <- dat[!is.na(dat$cr_intmin),]

# create chapter length variable
dat$rwe_chapter_length <- dat$page_right_end-dat$page_right_start
dat$lwe_chapter_length <- dat$page_left_end-dat$page_left_start

# ratio of chapter length
dat$rwe_lwe_chapter_ratio_log <- log(dat$rwe_chapter_length + 0.5) - log(dat$lwe_chapter_length + 0.5)

# remove those with NA in chapter ratio
dat <- dat[!is.na(dat$rwe_lwe_chapter_ratio_log),]

# logged crime ratio
dat$rwe_lwe_crime_ratio_log <- log(dat$extreme_crime_right + 0.5) - log(dat$extreme_crime_left + 0.5)

# get crime ratio on federal level as a separate variable and merge with the main data 
crime_fed <- dat[dat$jurisdiction=="bund",c("year", "extreme_crime_right", "extreme_crime_left")]
names(crime_fed) <- c("year", "extreme_crime_right_fed", "extreme_crime_left_fed")
crime_fed$rwe_lwe_crime_ratio_log_fed <- log(crime_fed$extreme_crime_right_fed + 0.5) - log(crime_fed$extreme_crime_left_fed + 0.5)
dat <- merge(dat, crime_fed, by="year", all.x=T)

# impute missing values in the logged crime ratio with value on the federal level
dat$rwe_lwe_crime_ratio_log_fedimp <- ifelse(is.na(dat$rwe_lwe_crime_ratio_log),
                                             dat$rwe_lwe_crime_ratio_log_fed,
                                             dat$rwe_lwe_crime_ratio_log)

# create bias measure (difference between logged chapter length ratio and looged crime ratio)
dat$bias_ratio_fedimp <- dat$rwe_lwe_chapter_ratio_log - dat$rwe_lwe_crime_ratio_log_fedimp

# create decade indicator
dat$decade <- NA
dat$decade[dat$year<1960]<-1950
dat$decade[dat$year>=1960 & dat$year<1970]<-1960
dat$decade[dat$year>=1970 & dat$year<1980]<-1970
dat$decade[dat$year>=1980 & dat$year<1990]<-1980
dat$decade[dat$year>=1990 & dat$year<2000]<-1990
dat$decade[dat$year>=2000 & dat$year<2010]<-2000
dat$decade[dat$year>=2010 & dat$year<2020]<-2010
dat$decade[dat$year>=2020 & dat$year<2030]<-2020

# merge with far-right polling data
polbar <- read.dta13("polbar_agg.dta")
polbar$jurisdiction <- tolower(polbar$jurisdiction)
dat$jurisdiction <- gsub("ü", "u", dat$jurisdiction)
dat <- merge(dat, polbar, by.x=c("jurisdiction", "year_pub"), by.y=c("jurisdiction", "year"), all.x=T)

## regressions models
# position
mod_1 <- lm(rwe_lwe_chapter_ratio_log ~ factor(cr_intmin)*polbar_fr_vote_pct, data=dat)
mod_2 <- lm(rwe_lwe_chapter_ratio_log ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction), data=dat)
mod_3 <- lm(rwe_lwe_chapter_ratio_log ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction) + factor(decade), data=dat)
mod_4 <- lm(rwe_lwe_chapter_ratio_log ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction) + factor(year), data=dat)
# bias
mod_5 <- lm(bias_ratio_fedimp ~ factor(cr_intmin)*polbar_fr_vote_pct, data=dat)
mod_6 <- lm(bias_ratio_fedimp ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction), data=dat)
mod_7 <- lm(bias_ratio_fedimp ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction) + factor(decade), data=dat)
mod_8 <- lm(bias_ratio_fedimp ~ factor(cr_intmin)*polbar_fr_vote_pct + factor(jurisdiction) + factor(year), data=dat)

## coefficient plot
ests_poll =  rbind.data.frame(avg_slopes(mod_1, variables="polbar_fr_vote_pct", by="cr_intmin")[,3:5],
                              avg_slopes(mod_4, variables="polbar_fr_vote_pct", by="cr_intmin")[,3:5],
                              avg_slopes(mod_5, variables="polbar_fr_vote_pct", by="cr_intmin")[,3:5],
                              avg_slopes(mod_8, variables="polbar_fr_vote_pct", by="cr_intmin")[,3:5])
ests_poll = as.data.frame(ests_poll)
names(ests_poll) = c("cr_intmin", "coef", "se")
ests_poll$cr_intmin = rep(c("Center-Left", "Center-Right"), 4 )
ests_poll$model = rep(rep(c("w/o Fixed Effects", "w/ Fixed Effects"), each=2), 2)
ests_poll$model = factor(ests_poll$model, levels=c("w/o Fixed Effects", "w/ Fixed Effects"))
ests_poll$outcome = rep(c("Position", "Bias"), each=4)
ests_poll$outcome = factor(ests_poll$outcome, levels = c("Position", "Bias"))

ggplot(data=ests_poll, aes(x=model, y=coef, group=factor(cr_intmin), shape=factor(cr_intmin))) + 
  geom_errorbar(aes(ymax=coef+1.96*se, ymin=coef-1.96*se), width=0, position = position_dodge(width = 0.5)) +
  geom_point(position=position_dodge(width=.5), size=2.5) + 
  facet_wrap(~outcome) +
  theme_bw() +
  geom_hline(yintercept = 0, linetype="dashed") +
  theme(legend.position = "bottom", text = element_text(size=18),
        panel.grid = element_blank(),
        axis.text.x = element_text(angle = 30, hjust = 1)) +
  xlab("") + ylab("Marginal Effect of\nFar-Right Polling") +
  scale_shape_manual(values=c(16, 17), name="")

ggsave("Fig4_5.pdf", width = 7, height=5)
